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Abstract 



Temperature stabilization of oil and gas wells is used to ensure stability 
and prevent deformation of a subgrade estuary zone. In this work, we 
consider the numerical simulation of thermal stabilization using vertical 
seasonal freezing columns. 

<Z2 A mathematical model of such problems is described by a time-dependent 

i ^**i temperature equation with phase transitions from water to ice. The re- 

sulting equation is a standard nonlinear parabolic equation. 
y—i Numerical implementation is based on the finite element method using 

^- the package Fenics. After standard purely implicit approximation in time 

l/~) and simple linearization, we obtain a system of linear algebraic equations. 

Because the size of freezing columns are substantially less than the size of 

^*-' the modeled area, we obtain mesh refinement near columns. Due to this, 

. we get a large system of equations which are solved using high performance 

computing systems. 

o 
m 
^ H 1 Mathematical model 

> 

We consider a mathematical model that describes the distribution of tempera- 

r> ture with phase transitions at a given temperature T* in a domain ft = £l~ UJ7 + . 

jrt Here £l + (t) is a domain with liquid phase, where the temperature is above the 

phase transition temperature 

n+(t) = {x\x e n, t(x, t)>r*} 

and £l~ (t) stands for a domain with solid phase, 

fi - (t) = {x\x e fl, T(x, t) < T*}. 

The phase transition occurs at a phase change boundary S — S(t). 

For simulation of heat transfer with phase transitions, the classical Stefan 
model is used [I], [2]. This model describes the thermal processes accompanied 



by a phase change media, absorption and release of latent heat. We use 

(a(cb)+p+Lcb') f~- ^gradp, g radTn-div(A(^)gradT)=0 ) (1) 

where L is thelatent heat of the phase transition, k stands for the absolute 
permeability of the porous medium, p, is the water viscosity and p is the pressure 
in the soil. 

For the coefficients of the equation, we have the following relations 

a(0) = P~ c ~ + <P(P +C+ - P~ c ~), 

and 

0, T<T*, 

1, T>T*, 

where p + , c + p~ , cr are the density and specific heat capacity of the melt and 
frozen zone, respectively. 

Since we consider the process of heat propagation in porous media, then for 
the coefficients we have: 

c~p~ = (1 - m)c sc p sc + mcipi, 

c + p + = (1 - m)c sc p sc + mc w p w , 

where m is the porosity. Indexes sc,w,i denote the skeleton of the porous 
medium, water and ice. For the coefficients of thermal conductivity in the melt 
and frozen zone, we have similar relationships 

A - = (1 - m)\ sc + m\i, 

A + = (1 — m)\ sc + m\ w . 

In practice, phase transformations do not occur instantaneously and can 
occur in a small temperature range [T* — A, T* + A] [T]. As the ^-function, we 
can take (f>A'- 



o, 




T < T* - A, 


T-T* 
2A 


+ A 


T* - A < T < T* + A, 


1, 




T > T* + A, 


o, 




T < T* - A, 


1 
2A' 


rp* 


- A < T < T* + A, 


o, 




T > T* + A. 



Then we obtain the following equation for the temperature in the domain $7: 
(a(<p A ) + PlLcf>' A ) f ^+MgradTj - div(A(<M)gradT) = 0. (2) 



The resulting equation pi) is a standard nonlinear parabolic equation. 

The equation (f2j) is suplemented with the initial and boundary conditions 



(*,0)=T , 


X £ 0, 


T = T C , 


x£T D 


on 


ier/] 



(3) 



Here Tjj is a place of contact with the freezing columns. 



2 Finite element realization 

The equation (pi) is approximated using a finite element method. We multiply 
the temperature equation by a test function v, and integrate it using the Green 
formula 

dT 

(«Oa) + Pi l 0'a)^7 v dx 

9t , (4) 

+ / (A(>a) gradT,grad«) da; = 0, W G H^(fi). 

Here iJ 1 (O) is a Sobolev space, which consists of the functions v such that v 2 
and \Vv\ 2 has a finite integral in the fl and Hq(Q.) = {v £ J? 1 (ft) : w|r D = 0}. 

To approximate in time, we apply the standard fully implicit scheme and 
use the simplest linearization (from the previous time level) 

(a(fi) + n L<t>t) vdx 

I <5) 

+ / (A(^)gradr n+I ,gradw)da; = 0. 

Jo. 

In such a way, we arrive at the following classical variational formulation of the 
problem: to find T £ H 1 ^), [T(x,t) - T c (x,t)] |r D G fl^(fi) such that 

f (a((f>l) + P iL^)T n+l v dx + f (A(^)gradT" +1 ,gradv)da; 

JO (6) 

{a(cj>l)+piLcf>'£)T n vdx, Vv£H 1 {n). 
The process of solving equation (TgJ) can be represented as follows. While 

"cur ^ t r> 



L cur ^ L max- 



2. Save a previous time lavel values T prev = T; 

3. Recalculate the ambient temperature T air : 

T a „ = 41sm((27r(W/86400 + 250))/365) - 10.2; 



4. If the temperature of the soil is less than T a i r , then the freezing columns 
turn on, else turn off; 

5. The temperature at the new time level are solved by a linear solver; 

6. Write the results to a file. 

Numerical implementation is performed using the Fenics package [5]. For 
results visualization, the values at each time level were recorded in vtk file 
format that was visualized using ParaView program. 

3 Numerical results 

As a model problem, we consider the process of thermal stabilization of the 
mouth of the oil or gas wells [3], [4]. The geometric domain was built using 
the Netgen mesher. The computational domain is shown in Fig. IT] and has a 
length of 40 meter in each direction, the well (radius 0.1 meter) is located in the 
middle of a field, where oil flows with a given positive temperature. A cement 
layer with the thickness of 0.2 meter is used for well heat insulation. 8 freezing 
columns with a radius of 0.05 meter are deepened to 14 meters around the well. 
The top sand layer has the thickness of 2 meters. Near the well, there is laid 
penopleks (10 by 10 m, thickness of 200 mm). 

Parameters of the problem are given in Table [TJ The computational grid 
contains 10,903,946 cells. 

The numerical results are presented in Fig. [3j Figure [5] illustrates the 
efficiency of the seasonal cooling devices, which accumulate the winter chill in the 
ground and provide additional bearing capacity in the summer. By numerical 
simulation of the temperature stabilization of soils using freezing system, we can 
conclude that the presence of freezing columns can reduce soil thawing around 
wells The calculations were performed using the NEFU computational cluster 
Ariane Kuzmin. The computation time using 64 processors (see Fig. p| was 
about 14 minutes, with 32 processors - about 21 minutes, and with 16 processors 
- about 40 minutes, which shows good efficiency of parallelization. 
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Figure 1: The computational domain 




Figure 2: Domain decomposition 



Notation Value Metrics 



Description 



J-cyl 


20.0 


degree 


To 


-5.0 


degree 


T, 


0.0 


degree 


L 


1.04e8 


J/kg 


CPsc 


2.17e6 


J/m3 


cpi 


2.42e6 


J/m3 


Cpsa 


1.34e6 


J/m3 


CPpe 


0.20e6 


J/m3 


CPce 


0.8e6 


J/m3 


^sc 


2.43 


W/(m degree) 


A; 


2.22 


w/(m degree) 


A sa 


0.47 


W/(m degree) 


Ape 


0.03 


W/(m degree) 


Ace 


0.21 


W/(m degree) 



Temperature of oil in the well 

Inital temperature 

Phase change temperature 

Latent heat of the phase transition 

Volumetric heat capacity of soil 

Volumetric heat capacity of water 

Volumetric heat capacity of sand 

Volumetric heat capacity of polystyrene 

Volumetric heat capacity of cement 

The thermal conductivity of soil 

The thermal conductivity of water 

The thermal conductivity of sand 

The thermal conductivity of polystyrene 

The thermal conductivity of cement 



Table 1: Problem parameters 




Figure 3: The temperature distribution after 5 years 
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Figure 4: Temperature after five years. Slice x — 20. Without freezing columns 
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Figure 5: Temperature after five years. Slice x = 20. With freezing columns 
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